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Abstract This paper describes a pilot study into the 
spectral energy distribution (SED) fitting and the 
derivation of physical parameters for 19 galaxies ob- 
served as part of the Great Observatories All-sky LIRG 
Survey (GOALS) survey as observed with the Spitzer 
Space Telescope. For this we have used the pan-spectral 
fitting tools developed in a series of papers by Dopita 
and his co-workers. We show that the standard Lee and 
Draine 'astronomical silicate' model cannot provide a 
good fit to the silicate absorption features as observed 
in the heavily dust-extinguished (Ay ~ 50mag.) star- 
bursts. We have derived an empirical fit to the 'star- 
burst silicate' absorption in these objects. This absorp- 
tion curve is consistent with the silicate grains being 
systematically larger in starburst environments than in 
the local Galactic interstellar medium. We demonstrate 
the sensitivity of the SED fitting to each of the fit- 
ted parameters, and derive these parameters for those 
galaxies which do not have an embedded AGN. This 
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technique is simple and provides reasonably robust and 
uniform parameters for the starburst, especially as far 
as the star formation rate, population of old stars, com- 
pactness of the starburst region and total foreground 
extinction are concerned. However, the chemical abun- 
dances and the optical extinction cannot be reliably 
determined by this analysis, and optical SEDs will also 
be required to provide a complete characterization of 
the starburst region and of the surrounding galaxy. 

Keywords ISM: dust -silicates — extinction, Flu - 
galaxies:general,star formation rates, spectral energy 
distributions, starburst-infrared:galaxies 



1 Introduction 

The spectral energy distribution (SED) of a galaxy is 
the observational quantity needed to derive most of the 
the fundamental parameters of that galaxy. From its 
study we can hope to determine the star formation rate, 
the quantity of older stars, and their star formation 
history, the amount of dust and its nature, and the 
chemical com position of both the stellar and inte rstellar 
components ( Dopita 2005a : Walcher et al. 2010h . 

For starburst galaxies, a wide variety of empirical, 
semi-emirical and purely theoretical techniques have 
been applied to the study of their SEDs, Fo r a re - 

" (l2010lh 



cent comprehensive review see IWalcher et al 



In the sem i-empirical mode l ling of Dale and hi s col- 
laborators (|Dale et al.ll200l[ bale fc Heloul 120021 ) . the 
SEDs of both disk and starburst galaxies were taken as 
a one-parameter family in terms of the dust tempera- 
ture, with more active star for mation being related to 
higher mean dust temperature. lLagache. Dole fc Puget 
( 2003!) derived an empirical sequence of SEDs in terms 
of the absolute luminosity of the galaxy. 
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Both of these semi-empirical approaches are valid 
to some extent, since IR luminous galaxies generally 
have greater rates of star formation, more compact star 
formation regions and higher luminosities than normal 
galaxies. Both of these models have found their greatest 
application in providing a simple one-parameter fit to 
the SEDs of normal disk galaxies. 

In starburst galaxies - to which we now confine 
our discussion - the geometry of the star forming re- 
gion is very important in determining the form of the 
IR SED at the wavelengths observed by the Spitzer 
Space Telescope. The SED in this wavelength region 
is dominated by the emission from the hotter dust 
grains located in and around the Hn regions associ- 
ated with the young stellar clusters, and by the poly- 
cyclic aromatic hydrocarbon (PAH) emission located 
in the dense photodissociation regions (PDRs) sur- 
rounding these Hll regions. PAHs are destroyed by 
ionizing photons, so cannot exist in the H II regions, 
but they require energetic UV photons absorbed in the 
2200A band (and shortwards) to excite the character- 
istic C-H, C-C and C- C-C skel e tal defo rmation modes 
of the PAH IR bands iDraine! (|2003al ). This spatial 
segregation of the PAH emission into the PDR is very 
clearly s een in the Spitze r Space Telescope image of 



M20 by iRho et al 



( 20061) . The spatial scale of the 
PDR regions from which most of the IR emission orig- 
inates is fairly small. For example, around an Hll re- 
gion with P/k ~ 3 x 10 6 cm~ 3 K, typ i cal of starbursts 
(|Kewley. Heisler. Dopita fe Lumsdenl l200lh . a PDR 



with Ay ~ 3 and T ~ 1000K (iTielens fc Hollenbach 
1985( 1 is only ~ 0.2pc thick! 

The star formation regions of starburst galaxies are 
both denser and of higher pressure than the more ex- 
tended gas and dust associated with the older stars. 
The truth of this can be sim ply demonstrated by con- 
sidering the lKennicuttl (|1998f l star formation law, which 
has been shown to apply both to normal disk galax- 
ies as well as starburst galaxies. This states that the 
local surface rate of star formation, Ssfr, is related 
to the surface density of gas, S gas , by a simple power 
law; Esfr °c ^ g & s - Thus, by definition, high local spe- 
cific star formation rates will always be associated with 
dense, high density regions with both high dust tem- 
peratures and high foreground extinction. However, it 
does not necessarily follow that the IR luminosity and 
the dust temperature are simply correlated, since the 
luminosity scales as the total star formation rate, while 
the the average dust temperature is determined by the 
average geometry of individual star formation regions. 
For example, a galaxy with many individual H II re- 
gions spread through a disk may have a high total star 
formation but such a galaxy would have a lower dust 



temperature than one in which the same star formation 
rate is crushed into a dense nuclear starburst. 

The simplest geome tries simply co-mix the star for- 
matio n and the dust ([Devriendt. Guiderdoni fc Sadat 
1999(1 . A somewhat more sophisticated a pproach was 



tak en by Ida Cunha. Chariot fc Elbad (|2008h . who used 
the ICharlot fc Falll (|2000h approach to dust absorption 
and emission. This recognises that the younger stars 
are more dust-embedded than the older stars, and that 
the 'birth cloud' emission is hotter than the more 'dif- 
fuse ISM' dust emission. However, the dust emission 
templates are rather simple, consisting of a PAH tem- 
plate and variable gray-body contributions. 

Spherical radiative transfer models have been much 
employed in the literature. This approach was pio- 



neered by Rowan -Robinson and his co-workers ( Rowan- Robinson] 
ll98Ct Il989l 1 19931) and developed to a high degree of so - 
phistication by Tagaki. Vansevicius fe Arimoto ( 20031) ; 
Tagaki. Arimoto fe Hanami (12003) and by Seib e nmor 



gen and his team(see ISiebenmorgen fe Kriigell ([20071 ) 
and references therein) . These have been applied to the 
determination of star formation rates in Starbursts. 

The difficulty with all of these models is that they do 
not contain significant Hll region physics. This leads 
to inconsistencies in the pressure and density of the 
gas, since the equilibrium temperature is not explic- 
itly solved for. Furthermore, they cannot reproduce 
the hot dust seen only in H 11 region. Thirdly, they 
fail to account for nebular line and continuum emis- 
sion processes. In starbursts, the Hn regions around 
clusters expand dynamically during the stellar lif etime . 
Efstathiou. Rowan- Robinson fe Seibenmorgen (2000) 



made models which explicitly followed the evolution 
of the giant molecular clouds owing to the ionization- 
induced expansion of the H II regions and the evolu- 
tion of the stellar population, following the Bruzual 
fc Chariot stellar popula tion synthesis models (see 
iBruzual fc Chariot ( 2003 ) for more up-to-date mod- 
els). In these models, the effect of transiently heated 
dust grains/PAHs on the radiative transfer, as well as 
multiple scattering, is taken into account. However, 
there is no gas physics. 

Dopita, Groves and colleagues have developed a ver- 
satile technique to model the pan-spectral energy dis- 
tributions of starburst galaxies. These methods cover a 
wide range of metallicities and physical parameters de- 
scribed briefly below. These models and their physical 
assu mptions have been fu l ly desc ri bed in a series of pa - 
pers (jPopita et al.ll2005bl l2006allbl ; iGroves et aill2007t) , 
and will be only brie fly recapitulated her e. The reader 
is encouraged to read lGroves et al] (|2007l ) in particular, 
since the models given there are used in throughout the 
remainder of this paper. 
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These models take into account all the atomic gas 
physics and radiative transfer through the Hn re- 
gions and their surrounding PDRs in the outward- 
only approximation. The effect of multiple scatter- 
ing is only included through approximate scattering- 
length methods within the Hn regions. The gas ion- 
ization state and temperature, and the grain tem- 
perature distribution and emission are explicitly com- 
puted at each radial step. All relevant emission lines 
and nebular continuum processes are all fully taken 
into account. The model dust has silicates, carbona- 
ceous grains as well as PAHs. The dust composition 
reflects the chemical composition of the gas. The 
amount of dust is calculated using solar depletion 
factors, and the size distribution of the grains is a 
power-law appropriate to a grain-shattering distribu- 
tion. The opti c al prop e rties of the grains are tak e n from 
Draine fe Led dl984l); IWeingartner fc Draind d200lh ; 
Li fc Drainel(|200ll) ; lDraine](|20Q3rj ) and lDrainel (|2003d) . 
The grain size spectrum is distributed into 80 bins 
for each grain type, and the emission spectrum of 
transiently heate d dust grains is comput e d foll owing 
the methods of Guhathakurta fe Drainel (1989) and 
Draine fc Lil (l200ll). For the PAHs, t he opacity is set 
according to the iLi fc Drainel ( 2001 ) approximation, 
and the ionization state and heating of the PAHs is ex- 
plicitly computed. The absorbed UV/visible radiation 
is re-radiated into a multi-component standard PAH 
emission band template in the mid-IR. 

Each H II region +PDR model provides the expected 
SED of an individual H n region at a chosen pressure, 
excited by a particular cluster, with a particular age. 
The art of our global SED modeling is to sum these 
models together so as to provide the self-consistent evo- 
lution of a whole ensemble of Hn regions of different 
ages and cluster masses within a dense high pressure 
region of the interstellar medium (ISM) . The dust tem- 
perature distribution which, for a particular choice of 
the metallicity, primarily determines the form of the IR 
SED, was shown to be set by a Compactness Param- 
eter. The concept of such a parameter was firs t de- 
veloped by Taeaki. Vansevicius fc Arimotol ( 2003 ) and 
Tagaki. Arimoto fc Hanamil (|2003h in the context of 
simple spherical radiative transfer models. In our case 
the compactness parameter C is defined in terms of the 
mean stellar cluster mass, M c \, and the pressure in the 
ISM, P Q , by: 
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In this paper we will apply these modeling techniques 
spectra drawn from a sub-set of the The Great Obser- 
vatories All-Sky LIRG Survey (GOALS) galaxies to in- 
vestigate the silicate absorption properties in starburst 
regions, and to self-consistently derive star formation 
rates and other physical parameters for those galaxies 
which do not contain an embedded AGN. 



2 Modelling the Starburst SED 



The SED fits devel oped in lDopita et all (|2005bl . 12006 
and 



Groves et al.l (|2007ft are based upon the stellar 



SEDs and the computed stellar mecha nical energy in- 



put f rom the STARBURST 99 code (|Leitherer et al 
19991) . We use a 1-D hydrodynamic code and our MAP- 



PINGS III code to calculate the emission line, PAH, 
dust a nd nebular continuum contributions. We also 
use the lFischera fc Dopital (|2003l I2004L 120051 ) turbulent 
screen models to compute the foreground attenuation 
at all wavelengths. 



The iGroves et al.l (|2007l) models used in this paper 
are based on the 2006 version o f the S TARBURST 99 
code. We have used a iKroupal (|2002l) broken power- 
law IMF between 0.1 and 120 M Q . Within the code, 
we use the standard combination of the Geneva and 



The complete set of parameters which enter into 
the SED fitting are described in the following section. 



Pado va stellar evolution tracks (|Vazquez fc Leitherer 
2005) to tabulate the mechanical energy delivered by 
the stars, both as stellar winds and as supernova explo- 
sions. In our SED fitting, the mechanical energy input 
is used to solve the 1-D hydrodynamical evolution of 
the H 11 regions, given a mean cluster mass and pres- 
sure in the ISM. Our MAPPINGS III code is then used 
to solve for the dust continuum, nebular continuum and 
emission line spectrum at each age, and the resulting 
spectra are then summed over all ages. To this H 11 re- 
gion spectrum we then add the spectra of ultra-compact 
H 11 regions associated with single stars, a dense circum- 
starburst molecular gas screen, the spectrum of older 
(10 < t < 100 Myr) stars which we assume to surround 
the starburst region, and fi nally a foreground t urbu- 
lent attenua ting dusty screen iFischera fc Dopital (|2003 . 
2004 120051 ) . These models are starburst-only in the 
sense that we do not include a stellar population more 
ancient than 100 Myr. 

The objective of this modeling was to provide a ge- 
ometry which is reasonably realistic, which correctly 
accounts for the evolution of the ensemble of H 11 re- 
gions, and which is characterized by as few physical 
variables as possible. Indeed, our final SED models are 
characterized by only seven parameters: 

• The current star formation rate, SFR. This is the 
total star formation rate in stellar clusters averaged 
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over the last lOMyr. This is a free parameter which 
determines the absolute scaling of the SED models. 
The metallicity of the ISM: Z/Zq. This is computed 
for the abundances of the STARBU RST 99 models: 
0.04, 0.1, 0.2, 0.4, 1.0 and 2.QZ/Z (jLeitherer et al 



1999). 



The compactness parameter C, defined above. This is 
computed in steps of 0.5 over the range 4 < logC < 8. 
The fraction, -Fuchii, of young (< l.OMyr) ultra- 
compact H II (UCHII) regions and their PDRs 



surrou nding individual massive stars (jDopita et al 



2006c). A value /uchii = 1-0 implies that within 



the first l.OMyr of evolution, the fraction of star for- 
mation rate producing UCHII regions is the same as 
the formation rate of stars in the massive clusters, i.e 
that 50% of all massive stars younger than 1.0 Myr 
are isolated and surrounded by ultra-compact H II re- 
gions. 

The optical extinction, Ay(SB)) associated with a 
dense turbulent foreground screen of dusty molecu- 
lar and H I gas surrounding the main star form ation 
complex (|Fischera fc Dopitall2003l I2004L l2005h . 
The fraction, F id of an older starburst stellar popu- 
lation with ages 10 < t < 100 My. This allows for the 
possibility that star formation has been ongoing for 
at least a dynamical timescale. A value of / id = 1 
corresponds to steady star formation, i.e. that the 
star formation rate in the period 10 < t < 100 Myr 
is the same as in the period of < t < 10 Myr. 
The optical extinction, Av( id)j o f a foreground tur- 
bulent attenuatin g dusty screen (jFischera fc Dopita 
20031 l2004l2005h . 



3 The GOALS Data 

GOALS combines data from the NASA Spitzer Space 
Telescope, Chandra X-Ray Observatory, Hubble Space 
Telescope (HST), and the Galaxy Evolution Explorer 
(GAL EX) observatories, together with ground-based 
data, into a comprehensive imaging and spectroscopic 
survey of over 200 low-redshift (z < 0.088), Luminous 
Infrared Galaxies (LIRGs)Q. 

These LIRGs represent a complete subset of the 
IRAS Revised Bright Galaxy Sample (RBGS), which 
comprises 629 extragalactic objects with 60/xm flux 
densities above 5.24 Jy, and Galactic latitudes above 
five degrees, and span the full range of galaxy types, 
AGN, and galaxy merger phases which are able to pro- 
duce objects so bright at f ar-IR wavelengths. The full 
GOALS sample is listed by lArmus et al.1 (|2009l ). where 



redshifts, luminosity distances and infrared luminosi- 
ties, Lir, are given using a cosmology with Hq = 70 
km s" 1 Mpc" 1 , il Y = 0.72, and ft M = 0.28. Through- 
out this paper we adopt the cosmological parameters 
H o =70.5 km s^Mpc" 1 , Q v =0.7 3, and n M =0.27 base d 
on the five-year WMA P results drlinshaw et al J 1 2 0091) . 

In a recent paper, Howell et al. ( 2010l ) has com- 
pared the UV and far-IR properties of these galax- 
ies. These authors find extreme specific star forma- 
tion rate (SSFR) in the GOALS, with a median value 
(3.9 x 10~ 10 yr" 1 ) - comparable to the highest SSFRs 
seen in the Spitzer Infrared Nearby Galaxies Survey 
sample. The GOALS LIRGs on average have greater 
IR excesses than would be expected based on their UV 
color, consistent with very large extinctions in the main 
star forming regions. Both of these properties are con- 
sistent with the GOALS sample representing the most 
extreme star formation environments in the local uni- 
verse. 



In another recent paper dU et al.llioioh have stud- 
ied the SEDs for 65 (U)LIRGs drawn from the GOALS 
sample, and spanning the entire electromagnetic spec- 
trum from X-ray to radio wavelengths. The ob j ective 
of this work is different from our study. U et al. ( 2010l ) 
used stellar population synthesis models to estimate 
the integrated stellar masses of these systems, and the 



Dale fc Heloul ([20021 ) templates to estimate the total 



1 see http://goals.ipac.caltech.edu 



IR luminosity, and to study a variety of color indices. 
Our objective in this paper is to test the starburst SED 
models, and to use these to derive what physical param- 
eters we can from the IR spectroscopy obtained with the 
Spitzer Space Observatory. 

The galaxies presented here are sources for which the 
full Spitzer Infrared Array Camera (IRAC) photometry 
and complete Infrar ed Spectrograph (IRS) spectroscopy 
(jHouck et all l2004h have been obtained. The IRS in- 
strument was used in both its high- and low- resolution 
modes. The wavelength bands covered by each of the 
four channels are 5.3-14.2 /im, 10.1-19.4 /xm, 14.3-36.3 
pLva, and 19-36.8 /im. The observed wavelengths were 
corrected to the wavelength in the rest-frame of the 
galaxy observed, and the fluxes were corrected to abso- 
lute vF v flux units (erg cm _2 s _1 ) using the luminosity 
distance. 

The data used here simply represents an exploratory 
subset of the GOALS galaxies, and is in no sense com- 
plete or even representative of the full GOALS dataset. 
We have simply selected the sources to be analyzed in 
this paper from amongst the 50 galaxies in the GOALS 
First Data Delivery. In order to ensure a large fraction 
of starburst-dominated sources, the galaxies were also 
selected to have strong PAH features. Double sources 
were deliberately excluded to avoid possible contamina- 
tion by a companion. Selection of sources also avoided 
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the obviously extended sources for the same reasons of 
beam-correction or confusion. In this sense, the galax- 
ies presented here are the morphologically more "sim- 
ple" sources in the sample, and are therefore presum- 
ably the most amenable to a simple integral SED anal- 
ysis. 

All the IRS spectra were extracted using the stan- 
dard pipeline, and flux calibrated assuming a point 
source. Thus these spectra represent "nuclear" spec- 
tra in all cases, and are not intended to provide an 
integrated MIR spectrum of the LIRG. In fact, many 
of the observed LIRGs a re somewhat resolved to IRS 
( Diaz-Santos et al.l I2009T ) . This means that since the 



4.1 Silicate Absorption 



slits are narrow, slit losses (and gains) are difficult to 
quantify. In some cases therefore, the orders do not 
match up at the boundaries of the spectra (at the 10% 
level) . For point sources this cannot happen, but for ex- 
tended sources it does because the Spitzer Space Tele- 
scope is diffraction limited and the spectrograph slits 
are narrow. The absolute flux calibration is based on 
point sources (stars), so that we automatically take 
into account slit and diffraction losses in the pipeline, 
but only for a point source. As a result, it was some- 
times necessary to normalize the spectra relative to each 
other. The low resolution data was helpful in this re- 
gard, since it overlaps with the higher resolution spec- 
tra. In the spectra presented here, we have taken the 
normalization wavelength to be 25/im, since this corre- 
sponds to the IRAS observation wavelength, thus facili- 
tating comparison with integrated fluxes measued with 
this instrument. Furthermore, the IRAC photometry 
had to be extracted and summed in a region which cor- 
responds as closely as possible to the IRS slit. Both of 
these steps have been taken for the data presented in 
this paper. In fitting the SEDs, both the high and the 
low resolution IRS spectroscopy is used. However, the 
resolution of the high resolution spectra is higher than 
the resolution of the models, and the observed line to 
continuum contrast is therefore higher than in our our 
model fits. 



4 Fitting the IR Spectra 

The IR fitting procedure is not a linear process, as each 
parameter affects a particular part of the IR spectrum, 
and in different ways. In what follows we will demon- 
strate this and show how most of the SED fitting pa- 
rameters can be derived in a fairly uncoupled way. How- 
ever before we discuss the detailed fitting, let us first 
discuss the issue that was encountered in first attempt- 
ing to fit the silicate absorption features at the nominal 
wavelengths of ~ 10 and 18/im. 



Our m odels - like most others- used the Li fe Drainel 
( 20011) "Astronomical Silicate" scattering and absorp- 
tion properties . The optical constants used in the 



Drainel (|2003bllch model do not depend on the labo- 
ratory measurements of any particular silicate material 
or any particular size distribution. Rather, they are 
empirically derived to fit the shape of the silicate ab- 
sorption seen the the Galactic neighborhood. In our 
fitting, we rapidly discovered that these parameters do 
not provide an adequate fit to either of the two absorp- 
tion bands in the most heavily reddened objects with 
Ay > 50mag. 

As an example, we show in panel (a) of Fig[T]a fit to 
ESO 507-G070 (IRAS 13001-2339). This fit implies a 
foreground extinction of Ay ~ 8mag. and a starburst 
extinction of Ay <~ 52mag. It is clear that the fit with 
"Astronomical Silicate" absorption fails in the following 
ways: 

• The peak of the 10/zm feature is too deep, and at the 
wrong wavelength. 

• The 18/im feature is not deep enough, is of the wrong 
shape, and its peak is also at the wrong wavelength. 

In order to better fit the silicate features, we also 
took an empirical approach, and adjusted the absorp- 
tion coefficients by hand until we were able to reproduce 
the depth, shape, and central wavelengths observed for 
silicate features in ESO 507-G070 (IRAS 13001-2339). 
In order to avoid changing the normalization, we kept 
the central absorption of the ~ 10/xm feature constrant. 
The resulting fit is shown in panel (b) of Fig [T] This 
implies a foreground reddening of Ay ~ 8mag. and a 
starburst extinction of Ay ~ 48mag. All spectral fea- 
tures in the spectrum are now well fit by this "Starburst 
Silicate" fit. 

The resulting absorption coefficients for our grain 
mixture are graphed in Figure [5] and are listed for gen- 
eral reference in Table Q] Note the shift of the peaks 
of the two absorption features to longer wavelengths 
(9.57/xm and 18.44/im respectively, instead of 9.11/im 
and 18.08/im ), the much stronger 18/xm absorption 
relative to the 10/zm feature, and the appearance of 
a "shoulder" on the absorption at 24//m. 

The form of the silicate absorption depends upon 
a number of variables. Both smooth structure of the 
observed silicate bands is an indication that the dust 
is highly amorphous rather than having a well-defined 
crystal structure. A consideration of the cosmic abun- 
dances and the condensation temperatures implies that 
the Mg and Fe silicates, such as amorphous olivines and 
pyroxenes is preferrred ( Ossenkopf. Henning fe Mathisl 
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Fig. 1 The comparision of the SED fit to ESO 507- 

G070 using (a) ^astronomical silicates" as defined by 

iDraine fc Led l) 19841 ) and (b) the modified "starburst sili- 
cate" absorption derived in this paper. The colored lines 
represent the IRS spectra, the blue points the IRAC pho- 
tometry points, and the stepped black line (histogram) is 
our SED model fit at each of the wavelengths computed in 
our model. 




log k (urn) 

Fig. 2 The comparison of the total (log 10 ) attenuation pro- 
duced by a grain m ixture including the "astronomical sili- 
cates" as defined bv IDraine fc Led (l984l l and one including 
the modified "starburst" silicate absorption derived in this 
paper. This was obtained by empirically fitting the observa- 
tions of the ULIRG galaxy ESO 507-G070 maintaining the 
normalization at the peak of thelO^tm absorption feature. 
Note that the "starburst silicate" curve shows a shift of the 
10/j.m feature to longer wavelength, a more pronounced and 
shifted 18/xm absorption feature, and a shoulder in the ab- 
sorption at about 2 4/xm. All of these are characteristic of 
larger Si grains (c.f. iDorschner et al l (|l995l ). figure 9). The 
absorption feature at log A ~ 0.78/mi is a water ice feature. 



19921) . The original IDraine fc Led (|1984l ) "astronom- 
ical silicate" absorption curve has a small 18/mi to 
10/xm absorption band ratio of 0.3 7. According to 



Ossenkopf. Henning fc Mathisl ( 19921) . values this low 



can be produced in pure Fosterite (Mg 2 Si04) or in ir- 
radiated olivine. By contrast, high ratios of the 18/im 
to 10/im absorption can be produced by different com- 
positions, such as a high mass fraction of Fe2SiC>4 or 
bronzite (Mgo.9Feo.iSi03. The absorption of bronzite 
(a form of py roxene) also displays the additional s houl- 
der at 24/xm (jOssenkopf. Henning fc Mathisl ll992). 
Using IR spectra of s ome of the ULIRGs ob served by 



Spoon et al.l (|2007n and llmanishi et al.l (|2007h with the 
Spitzer Space Telescope and performing radiative trans- 
fer usi n g the DUSTY c ode by Ivezic. Nenckova fc Elitzur^ 
<fl999lh ISirockv et al.l (j2008| l have used the silicate fea- 
tures to help constrain the dust chemistry. These 
observations also require dust with relatively high 
18/im to 10/im absorption ratio of the silicate fea- 
tures (around 0.5), and these authors fa vor cold 



dust model of lOssenkopf. Henning fc Mathisl ([19921 ) 
(rather than their warm dust m odel, or the "astro- 
nomical" silicate model of D raine ( Draine fc Lellll984l: 
IWeingartner fc Draine 2001 ). since thermal gradients 
in the dust layer can naturally produce the observed 
scatter in t he ratio of the 18/im to 10/zm absorptio n ra- 
tio and the lOssenkopf. Henning fc Mathisl (|1992l) cold 
dust model gives a high base value of this ratio. How- 
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ever, we have some problems with this interpretation. 
Whilst it may produce a change in the 18/im to 10/im 
absorption ratio between low-density ISM environments 
typical of disk galaxies, and the high density environ- 
ments characteristic of starbursts, it cannot at the same 
time explain the systematic shift to longer wavelengths 
of both of t he silicate absorption p rofiles which we have 



found here ( Dorschner et al. 19951) 



A simple explanation is that the grains in starburst 
environments are systematically larger than in the lo- 
cal ISM. Indeed, the rapid growth of grains in dense 
and cool molecular clouds is inevitable. The denser the 
environment, the more rapid will be the rate of grain 
growth, and we would expect that grains will grow to 
a size limit imposed by the depletion of the heavy ele- 
ments from the gas phase. If the molecular clouds are 
cool enough, grains can also grow by coagulation (al- 
though such grains will be very weak, and easily shat- 
tered if they find themselves in a warmer environment). 

The optical constants and the effect of size on olivine 
and pyroxene glas ses have been ca r efully measured in 



the laboratory by Dorschner et al. ( 19951 ). From this 



work it is found that for olivine glass particles, as 
their size increases, the 10/im feature initially shifts 
to longer wavelengths and then rapidly broadens when 
a > 1.25/im. At the same time, the ratio of the 18/xm to 
10/im absorption rapidly increases, the peak is shifted 
to longer wavelengths, and an absorp tion shoulder at 



24/xm becomes more pronounced see iDorschner et al 
(|l995h . figure 9). All of these properties correspond 
to the differences identified between the "astronomical 
silicate" absorption curve and our derived "starburst 
silicate" absorption. To provide a 18/im to 10/im ab- 
sorption band ratio of 0.63, and an absorption peak in 
the principal silicate feature at 9.57/im (as observed) 
would imply a mean grain size of 1.0 < a < 1.5/im 
in the starburst environment, a not unreasonable esti- 
mate. 

Thus, in conclusion, it seems almost certain that 
the change in the silicate absorption is the result of 
larger grain sizes in the starburst environment. How- 
ever, larger silicate grains are less effective in produc- 
ing optical extinction. Thus, for a given depth in the 
silicate features, the implied Ay will be less, but we 
would also expect that the implied column density of 
gas would be correspondingly larger. In what follows, 
these (uncertain) correction factors have not been taken 
into account in deriving the parameters of the starburst. 

4.2 Fitting the SED 

Fitting the SED to seven variables might be initially 
thought to be a very difficult procedure. Fortunately, 



however, most of these variables can be simply and un- 
ambiguously separated in these IR spectra, since they 
affect different spectral regions, and in different ways. 
In the following sub-sections we illustrate the sensitiv- 
ity of the model SEDs to each of these parameters in 
turn, by varying the parameters enough to make a sig- 
nificant change to the overall SED. The one parameter 
which cannot be determined at all from the IR SEDs 
is chemical abundance. In our fitting, we have varied 
the abundance in the range 0.4 to 2.0 solar, but this 
only significantly affects the scaling of the PAH emis- 
sion features. In principle, the abundance might be able 
to be estimated from the ratio of the PAH features to 
the continuum at 25 /imHowever, the uncertainties are 
at present too large to make this reliable. To derive 
accurate abundances will require optical SEDs and in 
particular, optical emission line fluxes, to be obtained 
for these galaxies. 

4.3 Star Formation Rate 

The fundamental and most important parameter to de- 
termine is the star formation rate (SFR). In the IR, the 
whole starburst region acts rather as a bolometer of the 
luminosity of young stars, since both the PAH and the 
dust continua scale as the star formation rate. Chang- 
ing the SFR simply scales the spectra up or down. This 
is shown in panel (a) of figure (3). This demonstrates 
that the star formation rate can be measured from these 
spectra (at virtually any wavelength) to an absolute ac- 
curacy of ±10%. 

We elect to determine the star formation rate by 
least-squares minimization of the difference between 
our model and the observations at all wavelengths ob- 
served by the IRS. The IRAC data points are not used, 
since these are also sensitive to the older stellar fraction, 
as discussed in the next sub-section. 

4.4 Old Stars 

The fraction of older stars, F \a, as mentioned above, 
is determined by the IRAC photometry. This works 
because the old stellar component only becomes sig- 
nificant in determining the SED at shorter wavelengths 
than the main PAH emission features. Panel (b) of Fig- 
ure (3) shows that the fraction F Q id is set by the nor- 
malization to the IRAC 3.5/xm and 4.5/im photometry, 
after the underlying star formation component has been 
set as described above. We always find an excess flux 
at the IRAC wavelengths that can be ascribed to the 
older stellar component, and we can fix this to an ac- 
curacy of about ±15%. For an unattenuated SED, the 
IRAC 3.6/zm to 4.5/im color is fixed essentially by the 



IRAS 05442+1732 

Sensitivity to SFR 
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Fig. 3 The sensitivity of the SED fitting to particular variables is demonstrated for one galaxy, IRAS 05442+1732. The 
stepped colored lines represent the models, and the black lines and points are the observational data. In panel (a) the star 
formation rate (SFR) is shown for three values, 27, 33 and 40 M0yr -1 . In panel (b) the old stellar population fraction, 
F id, is given fo the values 0.6, 0.9 and 1.2. In panel (c) we vary the compactness parameter logC from 5.5 through 6.0 to 
6.5. Finally, in panel (d) we adjust the fraction of ultra-compact HII regions, Fuchii, in the steps 0.0, 0.2 and 0.4. 
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slope of the Rayliegh- Jeans tail of the Black-Body dis- 
tribution. As we show below, a flattening of the IRAC 
3.6/imto 4.5/im color may be used to help constrain the 



foreground attenuation, A 



V(old) 



but this color index is 



not particularly sensitive as an attenuation measure. 

4.5 Compactness Parameter 

The compactness parameter, logC, a measure of the 
geometrical relation between the dust and the heating 
stars determines the temperature distribution of the 
dust. Higher compactness parameters move the peak 
of the far-IR emission to shorter wavelengths. How- 
ever, at wavelengths below about 30 /jm the slope of 
the dust continuum emission is essentially unchanged 
as the compactness parameter is increased. However, 
for compactness parameters above logC > 6 a flatten- 
ing of slope in the dust continuum starts to become 
evident. The main effect is to increase the 20-30/im 
continuum with respect to the PAH features, and to 
fill in the region between the principal PAH features at 
about 10 fim. This is shown in Figure (3) panel (c). 
We can probably determine logC to within ±0.3 from 
a combination of the ratio of the 7/im PAH features 
to the 9/Ltm continuum, the degree of veiling of the 10- 
14/Ltm PAH complex, and the slope of the 18-34/um dust 
continuum. 

4.6 Ultra-compact Hn Regions 

The Ultra-compact H n regions produce an SED which 
is strongly peaked near ~ 20/xm as a consequence of 
the very hot dust in the H II region itself, which strongly 
competes with hydrogen for the absorptio n of the ioniz 



ing ph otons coming from the central star (jDopita et al 
2006a ) . This ~ 20/im peak helps to distinguish UC 



H II regions from the the compactness parameter, since 
they change the slope in the 1.2 < log(A//xm) < 1.5 
wavelength region, but make little difference to the SED 
in the lO^m'TAH gap" (see Figure (3) panel (d)). We 
can determine -Fuchii to an accuracy of ±0.2. 

The presence of an active galactic nucleus (AGN) 
may sometimes be confused with a large fraction of UC 
H II regions, since both contain very hot dust giving 
a strong ~ 20/j,m peak. However, the presence of an 
AGN reveals itself in other ways, such as through a 
strong featureless continuum which washes out the PAH 
features, and through the high-excitation emission lines 
which appear in the IRS spectra. This issue will be 
further discussed below. 

4.7 Dust Attenuation 

Very high dust attenuation, equivalent to several mag- 
nitudes of optical extinction are required before ob- 
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Fig. 4 The sensitivity of the SED in fitting to the at- 
tenuation of the starburst, Av(sb) in the wavelength range 
covered by the IRAC and IRS observations. This is given 
for three diffferent cases, a low-extinction object, panel(a), 
an intermediate extinction case, panel (b) and a high- 
extinction case, panel (c). It is clear from the goodness 
of the fit in these examples that for AyfSB) < 3mag. the 
attenuation is essentially unmeasurable, while for v4 V (sb) > 
5mag. the attenuation can be estimated to an accuracy of 
~ ±5mag. 
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servable effects become visible at IRAC and IRS wave- 
lengths. Thus, the accuracy with which we can de- 
termine the foreground or the starburst reddening are 
dependent upon the absolute values of these quanti- 
ties. This is illustrated in Figure (4) in which three 
separate cases are shown. We conclude that when 
^4(V)sb < 3mag. the attenuation is essentially unmea- 
surable. When A(V)sb > 5mag the attenuation can be 
estimated to an accuracy of ±5mag. Note the flattening 
in the color index between he IRAC bands at 3.6/im and 
4.5/mi for the most heavily reddened values. This can 
be used to estimate the foreground reddening in those 
objects with A(v) \d > 5mag. However, the optical and 
UV SEDs are very much more sensitive to reddening, 
and these should be used in preference, if available. 



5 Results 

The results of our fitting for the 16 objects which could 
be well-described by the "pure starburst" modeling are 
shown in table [5] and figures [5] et seq. Note that all of 
these are rather good fits to the observed spectra, in- 
cluding the emission lines, although we would need to 
have chemical abundances derived from optical spectra 
of these objects in order to better quantify this state- 
ment. 

As might be expected for LIRGs, all of our galaxies 
turn out to have high compactness parameters (5.5 < 
logC < 6.5), indicating high pressures in the star for- 
mation region. From the definition of the compactness 
parameter given in Eqn. (1), if the mean cluster mass 
is logMci - 3 x 10 5 Af o , then 5.5 < log(P/ft) < 8.0 
(K cm -3 ). All of these objects are dense, high-pressure 
star formation regions with relatively "hot" dust. In- 
terestingly, the inferred star formation rate is not cor- 
related with the compactness of the starburst for this 
sample. 

We would not necessarily expect to find such a cor- 
relation, because compactness is description of the ge- 
ometrical relationship between the dust and the stars 
which heat it, while star formation rate is a global prop- 
erty of the galaxy. However, statistically for starburst 
galaxies a correlation between these qu antities should 
exist, and indeed, is hard- wired into the lDale fc Helou 
(j2002T ) template fitting. The physical reaso n for a cor- 
relation is the existence of Kennicutt (|l998h star forma- 
tion law. This shows that the surface rate of star forma- 
tion Esfr is related to the surface density of gas E gas 
through the relation Ssfr oc £ ga g. Thus object with 
high rates of star formation tend to be more compact, 
except in the cases wher e the star formation regio n is 
greatly extended. Indeed, iDiaz-Santos et al.l ([20091 ) has 



recently found that in many LIRGs the extended com- 
ponent of their MIR continuum emission originates in 
scales up to 10 kpc, and may contribute as much as 
the nuclear region to their total MIR luminosity. In 
addition, the more extended a galaxy is in the mid-IR 
the cooler far-IR colors it has. The absence of any cor- 
relation between star formation rate and compactness 
in our sample is probably more related to our selection 
for sources which are obviously compact, predicated by 
our need to start modeling objects which are morpho- 
logically simple. 

Although the objects show a great variety in the at- 
tenuation associated with the starburst (0 < A(v)sb < 
48 mag), the general extinction associated with the 
older population of stars (derived mainly from the 
color of the two IRAC photometry points, as described 
above) is much lower (0 < A(v) ga i < 8 mag). This 
is consistent with our expectation that regions of high 
specific star formation rate have high optical depths to 
dust absorption. 

The old stars fraction lies in the range 0.5 < F \a < 
3.5, consistent with starburst ages up to a few 10 8 yr, 
similar to the dynamical timescales of galaxies. The 
highest value is found in NGC 6701, which has low red- 
dening of both the starburst, and in the foreground. 

The fraction, Fuchii, of young (< l.OMyr) ultra- 
compact Hn regions is relatively low in most ob- 
jects ( 0.0 < Fuchii < 0.35), consistent with the 
expected lifetimes of the ultra-compact phase of evo- 
lution (< 1 — 2Myr). -Ftjchii reaches large values 
in only two objects, NGC 3221 (Ftjchii ~ 0.6) and 
NGC 5936(F UC hii ~ 3). It is possible that both of 
these objects may have a low-level AGN contribution, 
since AGNs also produce hot dust as described above, 
and it would therefore be difficult to distinguish from 
the contribution of UCHII regions in the ~ 20/im re- 
gion of the spectrum. NGC 3221 is nearly edge-on, 
so any AGN would be difficult to detect. However, a 
fairly bright and compact rad io source was observed by 
Irwin. Saika fc English] d2000l). and an X- r ay det ection 
was made by iRephaeli. Gruber fc Persia (|1995l ). sug- 



gesting that there may be a buried AGN in this ob- 
ject. However, NGC 5936 is a nearly face-on galaxy 
with a circum-nu clear stellar ring, and shows no sign 
of X-ray emission ( Rephaeli. Gruber fc Persiclll995l) . so 
this appears to be a bona-fide starburst with many sin- 
gle OB stars rather than having the star formation 
concentrated into clusters. -Fuchii ~ 3 implies that 
the total current star formation rate in this system is 
4 x 3.8 - WMoyr- 1 . 

In table [2] we have g i ven th e star formation rates in- 
ferred by iHowell et all (|2O10h for those galaxies which 
are in common. All our estimated star formation rates 
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Fig. 5 The fit to the observed SEDs of GOALS starburst galaxies. The colored lines represent the IRS spectra, the dark 
blue points the IRAC photometry points, and the stepped black line is our SED model fit given for each of the wavelengths 
computed in our model. Note that the issues associated with the normalization of the spectra affect the goodness of fit for 
some objects (for example, UGC3608 and UGC03351 on this panel, and NGC5936, IRAS F18292-3513 and ESO 286-G035 
on subsequent panels). 
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Fig. 6 As figure [5] 



13 




Fig. 7 As figure[5] 
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tend to be systematically lower, and for a few cases, the 
star formation rates we have derived are dramatically 
lower. The cause of this discrepancy is purely due to 
aperture effects. We have extracted spectra only from 
the I RS aperture, not the whole galaxy as is the case 
in the Howell et al. ( 2010l ) study. We must also bear in 
mind that the total instantaneous star formation rates 
for some galaxies are higher when we include the con- 
tribution of the UCHII regions. For example, including 
the UCHII regions gives a total star formation rate for 
NGC 5936 of 15M yr _1 , as noted above, much c loser to 
the 25.23M ? vr~ l inferred bv lHowell et al] (|2Q10h. Like- 
wise, the total star formation rate in NGC 3221 is 
increased to 6.1M yr _1 by the addition of the UCHII 
component. 

Finally, we want to draw attention to the somewhat 
peculiar object IRAS 06076-2139, shown in the last 
panel of Figure 8. Whilst the majority of the SED is 
well fit by our starburst model, there is a strong excess 
in the PAH emission at about 12.5/xm. This feature 
has been seen in spectral maps of NGC 7027, where 
it has b een interpreted as due to (uncharged) PAH 



clust e rs I Rapaciol i Joblin fc BoisselH2005t iBerne et al 



20081 lTielens| l2008|) . This type of PAH spectrum has 
been observed in other ULIRGS, for example IRAS 
13 416-6243. whic h has been dubbed a "class C" source 
bv lTielensI (|2008h . IRAS 06076-2 139 was the subject o f 
a detailed integral field study by lArribas et al.l (|2008). 
This system is a very close double system, with rel- 
ative velocities > 500km s , which have recently un- 
dergone a strong collision. This collision, and the resul- 
tant shocks may have something to do with the peculiar 
PAH spectrum. 

5.1 Objects with Embedded AGNs 

Of the 19 objects selected we found clear evidence for 
embedded AGN in only three objects. This is a surpris- 
ingly low fraction for such a luminous set of objects, but 
this simply reflects the pre-selection of sources towards 
those that contained strong PAH features which are 
likely to be modeled by simple starbursts. It will need 
a more complete and unbiased sample to determine a 
reliable starburst / AGN fraction. In our sample, the 
SED of MCG-03-34-064 is completely dominated by the 
bolometric luminosity of the AGN, and both the con- 
tinuum associated with hot dust and the strong high- 
excitation emission lines expected in an AGN are clearly 
seen, see Fig 8. 

Two other objects show a mixed starburst + AGN 
excitation, IRAS 08355-4944 and IRAS 23436+5257 
(see Fig 9). Both of these objects show an excess in 
the IRAC 4.5/xm band, and this excess extends all the 



way up to at least 20/im. In these objects the PAH band 
strengths are suppressed, presumably by the underly- 
ing hot dust continuum. There is also clear evidence of 
the high excitation lines expected in objects containing 
and AGN. We have estimated a star formation rate of 
- 85 Moyr^ 1 in IRAS 08355-4944, and - 50 Moyr" 1 
in IRAS 23436+5257, but both of these are overesti- 
mates of the true star formation rate, as the contribu- 
tion of the AGN to the bolometric luminosity cannot 
be properly disentangled in these objects. One possible 
way would to use only PAH fitting to estimate the un- 
derlying star formation, and then to add an AGN SED 
to the intrinsic starburst SED. This was not done here 
because we do not yet have good AGN te mplates. This 



approach is being currently developed by Pctric et al 
(|2010h . 



6 Conclusions 

In this paper we have shown that the IR SEDs of star- 
burst galaxies in the GOALS sample observed with the 
Spitzer Space Telescope can be modeled to a high degree 
of accuracy by the model s which we have developed in 
a series of recent p apers (|Dopita et al.ll2005bl [iooialibl ; 



Groves et al.l 12007). The observed SEDs can be used 



with these models to constrain the absolute star forma- 
tion rate, total extinction, and the ~ — 1, ~ — 10 
and ~ 10 — lOOMyr history of star formation. We can 
also determine a compactness parameter which is re- 
lated to the pressure in the interstellar medium around 
the starburst. All of these sources have high compact- 
ness, and star formation rates are high, ranging up to 
210 Moyr" 1 . 

We have also established that the silicate absorption 
in such dense starburst environments is fundamentally 
different than in the local ISM. Although this may be 
related to composition, we believe that it is more likely 
related to a population of grains with systematically 
larger sizes in the starburst environments. 
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Fig. 8 MCG-03-34-064, an AGN-dominated SED. The 
principal line identifications are marked, note that the slope 
of the two IRAC photometry points is reversed, silicate ab- 
sorption (if preset) is very weak, and that the PAH emission 
is practically swamped by the strong underlying continuum. 
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Fig. 9 The two examples of starburst + AGN SEDs found 
in our sample. 
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Table 1 Starburst Silicate Absorption compared with the Astronomical Silicate absorption in the 7.5-40 wav e length 
region. The column D&L84 is the attenuation of flux at wavelength A for Av = l.Omag using the iDraine fc Led i| 19841 ) 
'astronomical silicates'. The column labelled "Starburst Si" is the attenuation with our modified silicate absorption curve. 
The normalization of the two curves is identical at the peak of the 10/im feature. 
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Table 2 Fitted parameters for the 'pure' starburst galaxies. SFRj nt is the integrated star formation rate as derived for 
the galaxies in common with our sample bv lHowell et alj (|2010T ) 
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